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Preface 


The  purpose  of  this  research  was  to  investigate  and  compare  the  performance 
of  C-  and  O-grid  generation  methods  as  applied  to  predicting  the  flow  about  a  NACA 
0012  airfoil.  The  Beam  and  Warming  algorithm  was  employed  to  solve  the  Navier- 
Stokes  equations  as  applied  to  a  two-dimensional,  viscous,  compressible  flow.  Similar 
flow  conditions  were  establshed  utilizing  both  grid  generation  methodologies  and  the 
results  compared  to  experimental  data.  The  NACA  0012  airfoil  was  chosen  because 
it  has  been  extensively  studied  and  experimental  results  were  readily  available. 

Numerical  grid  generation  has  become  an  important  tool  used  in  the  numerical 
solution  of  partial  differential  equations.  The  general  method  involves  a  discretiza¬ 
tion  of  the  physical  domain  into  a  collection  of  grid  points.  The  differential  equations 
are  then  approximated  as  a  set  of  algebraic  equations  applied  in  this  domain.  In  sim¬ 
plest  terms,  the  grid  is  orthogonad  and  rectilinear.  In  the  case  of  a  curved  boundary, 
such  as  an  airfoil,  it  is  advantageous  to  consider  a  conformal  grid  which  greatly 
simplifies  the  application  of  the  boundary  conditions.  The  mapping  between  the 
physical  and  computational  domains  is  accomplished  by  means  of  the  transformar 
tion  metrics,  which  must  also  be  applied  to  the  field  equations  one  is  attempting  to 
solve.  The  accuracy  of  the  numerical  solution  is  therefore  influenced  by  the  type  of 
grid  chosen  and  the  associated  metrics.  Differing  methods  of  handling  the  boundary 
conditions  may  also  have  an  influence  on  the  numerical  solution.  For  this  reason, 
the  differences  in  the  metrics  and  the  application  of  boundary  conditions  are  also 
compared. 

I  would  like  to  thank  my  thesis  advisor,  Capt  Philip  Beran  for  his  continuing 
support  and  encouragment  during  this  research.  I  would  also  like  to  thank  Dr  Joseph 
Shang,  the  sponsor  of  this  research,  for  his  strong  support,  for  providing  computer 
resources,  and  most  especially  for  allowing  me  access  to  the  outstanding  group  of 
people  he  supervises.  Without  the  help  of  those  mentioned  this  thesis  could  not 
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have  been  accomplished.  Finally,  a  special  thanks  to  my  wife  and  children  for  their 
patience  and  moral  support  over  the  past  six  months. 
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Abstract 

The  purpose  of  this  study  was  to  investigate  and  compare  the  performance 
of  C-  and  O-grid  topologies  as  applied  to  numerically  predicting  the  flow  about 
a  NACA  0012  airfoil.  Both  types  of  grid  were  generated  using  a  hyperbolic  grid 
generation  code.  The  solution  of  the  flow  field  was  numerically  calculated  using  the 
Beam  and  Warming  approximate  factorization  method  to  solve  the  two-dimensional 
Navier-Stokes  equations  for  viscous,  compressible  flow. 

This  study  examined  the  performance  of  each  grid  for  a  range  of  Mach  number 
(0.3  to  1.1)  and  angle  of  attack  (2  to  12  degrees).  The  Reynolds  number  was  approx¬ 
imately  2,000,000,  with  variations  to  exactly  match  experimental  flow  conditions. 
The  numerical  inquiry  consisted  of  two  parts.  The  first  part  was  a  grid  sensitivity 
study  using  the  C-grid.  Progressively  finer  grids  were  employed  until  the  respective 
numerical  solutions  converged.  The  second  part  was  a  parametric  study  using  the 
O-grid.  The  trailing  edge  curvature  was  reduced  from  0.005c  to  0.0015c  and  the 
numerical  results  compared  with  those  of  the  C-grid. 

The  numerical  solutions  obtained  with  both  the  C-  and  O-grid  configurations 
compare  favorably  to  both  experiment  and  previous  numerical  results.  The  numerical 
results  produced  indicate  that  both  grids  are  satisfactory  in  a  general  application. 
However,  specific  weaknesses  arise  for  each  type  of  grid  in  certain  situations.  In  these 
instances,  there  may  be  an  advantage  in  selecting  a  particular  type  of  grid. 

In  the  subsonic  range,  the  O-grid  is  superior  to  the  C-grid  in  determining 
the  pressure  coefficient  in  the  vicinity  of  the  leading  edge.  The  O-grid  displays  an 
advantage  when  large  flow  gradients  appear  in  the  vicinity  of  the  trailing  edge,  such 
as  the  formation  of  a  lambda  shock.  The  O-grid  also  attains  convergence  for  certain 
critical  cases  (high  Mach  number,  Reynolds  number,  or  angle  of  attack)  for  which 
the  C-grid  solution  diverges. 


The  C-grid  displayed  several  advantages  over  the  O-grid.  One  advantage  is 
ease  of  application.  It  is  generally  easier  to  construct  a  C-  grid  to  a  given  set  of 
specifications  than  an  O-grid.  The  C-grid  handles  the  flow  resolution  better  in  the 
wake,  especially  when  the  flow  gradients  are  concentrated  directly  behind  the  airfoil. 
In  this  region  the  C-grid  displays  a  considerably  higher  concentration  of  grid  points. 
The  C-grid  also  predicts  the  pressure  distribution  at  transonic  speeds  and  the  mid- 
chord  shock  location  more  accurately  the  O-grid. 


COMPARISON  OF  C-  AND  O-GRID  GENERATION  METHODS 
USING  A  NACA  0012  AIRFOIL 


I.  Introduction 

The  purpose  of  this  study  is  to  investigate  and  compare  the  performance  of 
two  differing  grid  generation  methods  with  varying  topologies  as  applied  to  the  flow 
about  a  NACA  0012  airfoil.  The  numerical  method  of  solution  is  based  on  the  Beam 
and  Warming  approximate-factorization  algorithm  applied  to  the  two-dimensional, 
mass-averaged,  Navier-Stokes  equations  for  compressible,  viscous  flows.  While  this 
algorithm  is  applicable  to  both  steady  and  unsteady  flows,  the  solutions  considered 
in  this  study  are  all  steady  state. 

Numerical  grid  generation  has  become  an  important  tool  used  in  the  numeri¬ 
cal  solution  of  partial  differential  equations  (PDBs).  The  general  method  involves  a 
discretization  of  the  physical  domain  into  a  collection  of  grid  points.  The  differential 
equations  are  then  approximated  as  a  set  of  non-linear  algebraic  equations  applied  to 
this  domain.  In  simplest  form,  the  grid  is  orthogonal  and  rectilinear.  In  the  case  of 
a  curved  boundary,  the  difficulty  of  extrapolating  data  at  gridpoints  adjacent  to  the 
boundaries  is  introduced.  Not  only  is  this  extrapolation  cumbersome  and  difficult, 
it  also  imposes  a  reduction  in  the  order  of  accuracy  for  the  boundary  conditions, 
thus  degrading  the  fidelity  of  the  solution  throughout  the  field.  To  overcome  these 
difficulties  for  the  cue  of  a  curved  or  irregular  boundary,  such  as  am  airfoil,  it  is 
advantageous  to  consider  a  conformal  grid  that  simultaneously  eases  the  implemen¬ 
tation  of  the  boundary  conditions  amd  avoids  the  loss  of  accuracy  associated  with 
the  aforementioned  method. 


Grid  generation  methods  can  be  classified  into  three  major  catagories  (6:521): 
complex  variable  methods,  algebraic  methods,  and  differential  equation  methods. 
The  differential  equation  method  is  the  most  common  and  powerful,  especially  for 
complex  geometries.  These  methods  generally  employ  solutions  of  Laplace’s  or  Pois¬ 
son’s  equations  to  specify  the  coordinate  transformation  between  the  physical  and 
computational  domains.  The  PDGs  chosen  are  generally  elliptic  or  hyperbolic,  al¬ 
though  some  recent  investigations  have  explored  parabolic  grid  generators.  The 
elliptical  methods  have  a  strong  applicability  to  internal  flows  and  also  to  external 
flow 8  where  multiple  bodies  are  present.  Hyperbolic  methods  are  generally  easier  to 
apply  and  are  applicable  for  external  flows  where  there  is  no  constraint  of  the  outer 
computational  domain.  The  grid  generation  code  employed  in  this  study  was  written 
by  Barth  and  Kinsey  (9)  and  is  a  modification  of  the  hyperbolic  method  described 
by  Steger  and  Chaussee  (15). 

Two  different  types  of  grid  are  considered  in  this  work:  a  C-grid  and  an  O-grid. 
The  C-grid  envelopes  the  airfoil  in  C-shaped  loops  beginning  and  terminating  in  the 
far  wake.  The  branch  cut,  or  the  line  along  which  the  numbering  of  the  grid  lines 
begins  and  ends,  therefore  exists  in  the  wake  and  extends  from  the  trailing  edge  to 
the  far  field  boundary.  The  O-grid  has  oval-like  loops  which  completely  encircle  the 
airfoil.  The  cut  line  for  the  O-grid  is  arbitrary,  and  for  this  study  extends  from  the 
far-field  boundary  in  front  of  the  airfoil  to  the  leading  edge. 

Of  particular  interest  when  comparing  two  different  grids  for  a  given  applica¬ 
tion,  aside  from  rectifiable  differences  such  as  grid  spacing  and  resolution,  are  the 
transformation  metrics,  which  provide  the  mathematical  map  from  the  physical  to 
the  computational  domains.  Since  the  metrics  appear  in  the  differencing  equations 
to  be  solved,  they  have  an  impact  on  the  local  truncation  error.  Also  of  interest  are 
any  differences  due  to  the  varying  means  of  applying  the  boundary  conditions. 

This  study  consists  of  two  basic  parts.  The  first  consists  of  a  grid  refinement 
study  using  the  C-grid.  Progressively  finer  C-grids  were  employed  for  selected  con- 


ditions  (Mach  number  ranging  from  0.3  to  0.78;  angle  of  attack  between  2  and  4 
degrees)  until  the  solutions  became  invariant.  To  incorporate  the  O-grid  generator 
for  the  NACA  0012  airfoil,  it  is  necessary  to  truncate  the  trailing  edge  and  fit  it 
with  a  circular  arc.  As  the  truncation  point  is  moved  aft,  the  resulting  trailing-edge 
radius  of  curvature  (TBC)  becomes  smaller.  The  second  part  of  the  study  involved  a 
parametric  study  based  on  the  TEC.  Values  of  TEC  ranging  from  0.005c  to  0.0015c 
were  considered.  The  results  are  compared  with  those  obtained  using  the  C-grid. 

Data  analyzed  for  comparative  purposes  included  the  lift  coefficient,  C/,  drag 
coefficient,  Cdl  and  pressure  coefficient,  0p.  To  examine  convergence  rate  differences, 
Ci  and  Cd  are  plotted  versus  iteration  number.  Contours  of  Mach  number,  pressure 
coefficient,  and  vorticity  are  presented.  Contours  of  the  Jacobian,  aspect  ratio, 
and  orthogonality  are  presented  and  examined  to  determine  structural  differences 
between  the  two  grid  types. 

Both  versions  of  the  Navier-Stokes  code  (for  C-  and  O-grids)  were  authored  by 
Dr  Miguel  Visbal  of  the  Wright  Research  and  Development  Center  (WRDC)  Flight 
Dynamics  Lab  (FDL)  (16)  (17).  The  code  is  an  implementation  of  the  implicit, 
approximate-factorization  algorithm  developed  by  Beam  and  Warming  (4).  The 
turbulence  model  is  a  modification  of  the  algebraic  eddy  viscosity  model  developed 
by  Baldwin  and  Lomax  (3)  (See  AppendixJ)). 

Computer  resources  required  were  provided  by  FDL  and  included  accounts  on 
the  Cray-Vax  and  the  Cray  XMP.  Data  plots  were  performed  on  the  FDL  Prime 
computer.  Selected  data  runs  were  performed  on  the  AFIT  Stellar  computer  for  the 
dual  purpose  of  comparison  and  code  checkout. 


II.  Analysis 


2. 1  Governing  Equations 

The  fundamental  equations  of  fluid  dynamics  are  based  on  the  principles  of 
classical  mechanics  and  thermodynamics:  the  conservation  of  mass,  momentum,  and 
energy.  These  equations  Me  often  grouped  together  and  referred  to  as  the  Navier* 
Stokes  equations.  The  equations  are  presented  in  an  Bulerian  formulation  using 
vector  notation,  and  are  valid  for  unsteady,  compressible,  and  viscous  flows  (13::11): 
Conservation  of  Mass  (Continuity) 

^  +  y(,V)  =  0,  (1) 

Conservation  of  Momentum  (Newton’s  Second  Law) 

+  y .  =  *>/+ v  z,  (2) 

Conservation  of  Energy  (First  Law  of  Thermodynamics) 

OEt  *  ■*  dQ  #  ^  ^  «*  # 

+  (v  ■  V)B,  = -£ + p}  y  +  v.(e  vo-v  i  (3) 


The  stress  tensor  <l  is  expressed  in  terms  of  the  fluid  pressure  and  the  viscous  stress 
tensor,  r. 


o  a  —pi  +  r.  (4) 

Other  variables  ore  defined  in  the  list  of  symbols  at  the  beginning  of  this 
report.Two  additional  relations  are  required  to  close  the  equations.  These  are  the 


equation  of  state  for  a  perfect  fluid 


p  =  pHT  s  pe(7  -  1), 


and  Fourier’s  law  of  heat  conduction  that  relates  the  heat  flux  vector  to  the 


temperature  gradient 


q=  -kVT, 


where  the  thermal  conductivity,  k,  is  a  function  of  the  temperature. 


If  a  Newtonian  fluid  is  assumed,  then  Stokes’  law  relates  the  viscous  stresses  to 


the  rate  of  strain 


r  =  A( V  •  V)I  +  p(VV  +  ( Wf), 


where  Stokes’  hypothesis  3A  +  2fi  —  0  essentially  equates  the  mean  pressure  in  a 
deforming  viscous  fluid  and  the  thermodynamic  pressure  (20:70).  The  total  internal 
energy,  Et,  can  be  related  to  the  specific  internal  energy,  e,  and  the  velocity  (in  two 
dimensions): 

(8) 


To  facilitate  the  solution  of  the  NS  equations,  they  are  nondimensionalized  (see 
Appendix  A)  and  are  given  below  in  conservation-law  form  for  a  two-dimensional, 
rectilinear  coordinate  system.  The  interned  heat  generation,  Q,  is  assumed  constant 


in  time. 


dp  d(pu)  d(pv) 
dt  *  8x  *  dy 

d(pv)  .  d(pua  +  p-r„)  ,  d(puv-Tty)  T 

sr  +  — Tx — +  — ai - pU 
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(11) 


d(pv)  ,  d{pv*  +p-T¥y)  t  d(fmv  -  tvi)  _  - 

dt  +  ^  +  di 


dEt 

dt 


+ 


-b + jr  }  +  a(a«tg-»/«.-y.!g±») 

OX 

djEjV  +  pv  -  HTSV  -  VT„  +  fr)  _ 


(12) 


The  components  of  the  non-dimensional  stress  tensor  fire  given  by 


+  6,^ 


dtt* 

dik 


(13) 


These  equations  are  applicable  to  sin  orthogonal,  rectilinear  coordinate  system.  To 
transform  the  equations  so  that  they  apply  to  a  general  curvilinear  coordinate  system, 
the  chain  rule  is  employed.  The  transformation  is  from  the  physical  domain  (x,y) 
to  the  computational  domain  (£,17),  where  the  partial  derivatives  take  the  following 
form  (13:252) 

t  =  f(*»y)  (14) 


V-V(xty). 


(15) 


The  direct  application  of  the  chain  rule  then  produces  the  derivative  operators 


d_ 

dx 


d_ 

dy 


d(  d  dr)  d 
dx  d(  +  dxdr\ 
,  d  d 

t'dt+^dr) 
d£  d  di)  d 
dyd(  +  dydi) 
,  d  d 

+  I,l'at) 


(16) 


(17) 

(18) 
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The  metrics  are  determined  in  the  following  manner 


di  =  tsdx  +  ivdy 

(19) 

dr)  =  r^dx  +  r)ydy. 

(20) 

In  matrix  form 


and 


di' 

dx 

dr) 

n* 

Vy 

.  ^  . 

dx 

*( 

di' 

dy 

4 

.  * 

dr>. 

-i 

-| 

-1 

r  i 

is 

iy 

*•» 

=  7 

yrt 

ns 

n, . 

.  * 

Vn  . 

-ye 

*«  . 

Therefore,  by  solving  the  transformation  identity  (19:141) 


6 

-l 

1  0 

»?* 

n, . 

i— _ 

0  1 

•  ■ 

and  defining  the  Jacobian  of  the  transformation 


j  W,V) 

d(*,y) 


s  i 

d*  a* 


—  —  £,0* 


(21) 


(22) 


(23) 


(24) 


(25) 
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J  - 


9s  9s 

Sf  Sf 


(26) 


the  following  relationship  between  the  metrics  can  be  established: 


6  =  VnJ 

ns  =  -y(J 


(y  —  —XffJ 

n>  =  -x(J. 


(27) 

(28) 


The  NS  equations  may  then  be  written  in  terms  of  the  general  curvilinear  coordinates 
(18:5) 

dU  ,  dF  ,  dG  dF  dG  n 

at  +  *°»  t29) 

where  body  forces  are  omitted  and 


U  =  (p,pv,pv,pey 

l 


F  = 


G  = 


pu 

pv?  -  Ttt 
puv  —  Tsy 
(pe  +p)u-Fi 

pv 

PUV  -  TZy 

pv *  -  Tn 
^  (pe  +  p)v  -  Gt  ) 


(30) 


(31) 


(32) 
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To  incorporate  turbulence,  the  viscosity  is  separated  into  the  molecular  viscos¬ 
ity,  /*,  and  the  turbulent  eddy  viscosity,  e.  Then  Stokes’  hypothesis  is  extended  to 
include  the  turbulent  eddy  viscosity  as 


A<  =  -|(M  +  0-  (33) 

The  molecular  viscosity  is  calculated  from  Sutherland’s  formula.  The  turbulent  eddy 
viscosity  is  obtained  via  the  eddy  viscosity  model  of  Baldwin  and  Lomax  (3)  (See 
Appendix  D).  In  a  similar  fashion,  the  thermal  conductivity  is  separated  into  a  nom¬ 
inal  component,  k,  and  a  turbulent  component,  &t,  by  defining  the  nominal  and 
turbulent  Prandtl  numbers  respectively 

k'=c*-k-  (34) 

The  viscous  stress  and  heat  conduction  terms  can  then  be  written  in  the  following 
form  (18:6) 


= 

-3\us  +  A,(«,  +  vv) 

(35) 

3= 

+  vy) 

(36) 

ryy 

= 

+  A,(u,  +  v„) 

(37) 

Pi i 

= 

UTSS  +  t>rxy  -  qs 

(38) 

Gt 

SB 

urxy  +  vrn  -  qv 

(39) 

P 

S 

p(lf  -  l)(e  -  ^(u7  +  v2)) 

(40) 

= 

(«) 

<h 

= 

(42) 

The  molecular  Prandtl  number,  Pr,  the  turbulent  Prandtl  number,  Prit  and  the 
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specific  heat  at  constant  pressure,  Cp,  are  taken  as  constants: 


Pr  = 

a  .72 
k 

(43) 

Pft  = 

£Cp  e *  q 
k 

(44) 

°p  ~ 

a  1,000  Joules/Kg-K. 

7-1 

(45) 

The  NS  equations  can  be  cast 

in  strong  conservative  form  as  follows  (18:7) 

dtJ  dF  9G_ 
dt  +  d(  +  dr) 

(46) 

where 

«  =  vi 

(47) 

f.  U.F  +  (,G) 

^  J 

(48) 

A  (TtsF  +  tyG) 

G  -  , 

(49) 

To  facilitate  the  numerical  implementation,  the  strong-conservation  form  of  the  equa, 
tions  is  rewritten  as  (18:8) 
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dU  dEx  dE2 
—  +  — -  +  — - 
dt  di  dv 


dVi(uA)  |  dHVA)  ,  ,  dW2(U,Un) 

r\  +  I  ^  . 


(50) 


di 


d( 


di) 


dr) 


where 


Ex 


\ 


i  pvii  +  isV 

3  pvli  +  (yP 

K  (P  +  Pep  > 


/ 


pV 


\ 


P«v  +  r)sp 

pt/V  +  rjyP 

k  (P  +  P«)^  > 


/ 


0 


\ 


&i«e  +  6jvf 
i»at4f  +  b3V( 

bivtit  +  &a( t/«(  +  «V()  +  +  btT{  t 


/ 


0 


\ 


Citt,,  +  c2vn 

C3«7  +  CiVr, 

k  C\UVf)  +  CjUVij  +  +  c4t/t/,,  +  caTa  t 


Ci«c  +  C3t>{ 
CjUf  +  C4V( 


\  Ci«u^  +  c3vuf  +  C3UV4  4-  ctW(  +  c&T(  I 


(51) 


(52) 


(S3) 


(54) 


(55) 
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0 


IV,  =  I 

dit»„  +  djt»„ 

J  ./ 

djtt,  +  d3t/„ 

k  +  da(vti,,  +  at/,,)  +  d3vvn  +  diTn.  j 

U  and  V  denote  the  contiavariant  velocities: 


(56) 


U  =  tsu  +  t,v  (57) 

V  =  ffett  +  V'.  (58) 

The  viscous  coefficients  are: 

*  =  -|^(|e+«)  (59) 

(60) 

63  =  (61) 

=  *(£  +  £)@+S)  (62) 

Cl  =  (3^* +^,v)  (63) 

=  |^»  (|C.*V  “  (64) 

°3  =  +  (65) 

c<  =  2*1  (6*fc  +  (66) 

c*  =  — (p^  +  jp^r)  (^*^  4“  ^ty)  (67) 

*  =  -!*«(!<+<)  w 

da  =  -^A  fTjst),  (69) 

*  =  (™) 

•  *  “  ■>  (£  +  ^7)  +  tf)  ■  <71> 
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The  NS  equations  written  in  the  above  form  are  valid  for  unsteady,  viscous 
flows.  The  governing  equations  form  a  hybrid  hyperbolic-parabolic  system.  Their 
parabolic  nature  arises  from  the  second-order  derivatives  in  the  momentum  and 
energy  equations.  The  continuity  equation  contains  only  first-order  terms  and  is 
hyperbolic  (5:27).These  characteristics  of  the  system  of  partial  differential  equations 
dictate  the  manner  in  which  the  boundary  conditions  may  be  applied  (13:312).  See 
Appendix  C  for  an  explanation  of  the  boundary  conditions. 

2.2  Grid  Analysis 

As  previously  mentioned,  the  grids  used  for  this  study  were  of  two  types:  C-grid 
and  O-grid.  Both  were  numerically  generated  using  the  hyperbolic  grid  generation 
method  of  Steger  and  Chaussee  (15)  as  incorporated  in  the  code  written  by  Barth 
and  Kinsey  (9).  An  advantage  of  employing  a  hyperbolic  grid  generator  is  that  only 
the  inner  boundary  need  be  specified  (6:530),  i.e.,  a  point  distribution  is  specified 
on  the  airfoil  surface  only.  The  solution  then  proceeds  outward  from  the  surface  in 
a  time-like  manner. 

The  algorithm  developed  by  Steger  and  Chaussee  utilizes  the  definitions  of  the 
transformation  Jacobian  and  cell  orthogonality  to  devise  a  mapping  scheme  from  the 
physical  to  the  computational  domain  (9:3): 


+  VM  =  0 

Orthogonality 

(72) 

3  =  3^}  “  “  *h(s 

Jacobian 

(73) 

3 1  *  tjf# - w 

Inverse  Jacobian. 

(74) 

Here  the  Jacobian  represents  the  cell  area  ratio  between  the  computational  and 


physical  domains,  and  is  variable  throughout  the  mesh.  A  typical  transformation  is 
illustrated  in  Figure  1.  The  area  integrands  in  the  the  physical  and  computational 
domains  can  be  written  as  follows  (9:3): 


dA  =  dxdy  =  J~ld(dr} 

(75) 

Jdxdy  =  d£dr). 

(76) 

Numerically,  A(  and  At)  are  arbitrary  and  are  most  conveniently  choBen  to 
equal  one.  In  this  case,  the  inverse  Jacobian  approximates  the  physical  cell  area. 
The  orthogonality  constraint  and  inverse  Jacobian  are  locally  linearized  using 

x  =  x°  +  Ax  (77) 

y  =5  y°  +  Ay,  (78) 


where  x°  and  y°  denote  a  known,  near-solution  state.  The  linearization  then  proceeds 
as  follows  (15:433): 

x^xn  =  (x°  +  Ax)((x°  +  Ax\ 

=  x°(x°  +  x°(Axn  +  x°Ax(  +  0(  A3) 

=  f(*-*°)v +<(*-*’)< 

=  x\xn  +  -  x\x°n  (79) 
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l»)  Physics!  P!*ne 


Figure  1.  Transformation  from  Physical  to  Computational  Plane  (Ref:  IS) 
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The  resulting  locally  linearized  system  is  then  (9:4) 


i««  -* 


or  alternately 


x 

y 


* 


x 


x 

y 


o 

v  +  v° 


(80) 


AR<  +  BIir,=  F. 


(81) 


The  finite-difference  scheme  is  central  differenced  in  £  and  forward  differenced 
in  f).  Typically,  x|  and  are  approximated  as  follows  (9:6) 


x°t  = 


= 


2A£ 

1/t+ij  ~  ife-u 
2A£ 


(82) 

(83) 


while  s"  and  y?  are  extracted  from  the  orthogonality  constraint  and  the  Jacobian. 


*;  = 


*s  = 


(if)’  +  (x\? 

x°(V° 

wr+W 


(84) 

(85) 


There  are  certain  geometric  conditions  for  which  the  hyperbolic  algorithm  as 
described  has  difficulties.  These  are  pronounced  slope  discontinuities  and  regions  of 
reverse  or  concave  curvature.  These  problems  may  produce  unwanted  grid  charac¬ 
teristics  if  the  discontinuities  or  regions  of  reverse  curvature  propogate  into  the  grid 
interior.  By  adding  explicit  and  implicit  dissipation  these  difficulties  can  be  partially 
circumvented.  The  resulting  hyperbolic  equations  are  utilized  in  the  numerical  code 
of  Barth  and  Kinsey.  The  reader  is  referred  to  (9)  for  a  development  of  the  delta 
formulation  of  the  algorithm  employed. 
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Since  the  local  truncation  error  will  in  part  depend  upon  the  discretization  of 
the  coordinate  system  (i.e.,  the  transformation  metrics),  it  is  useful  to  have  an  alter¬ 
nate  means  to  examine  a  grid  in  terms  of  these  parameters.  A  physically  meaningful 
reparameterization  of  the  Jacobian  matrix  is  presented  by  Kerlick  and  Klopfer  (8). 
The  metric  tensor  is  introduced  as  an  intrinsic  property  of  the  transformed  geometry. 
The  metric  tensor  is  defined  as  (8:788) 


dxk  dxl 

a”  =  <,l3pSfJ' 


(86) 


where  6 u  represents  the  Kronecker  delta  function 


Ski  — 


0 

1 


ilk  £  l 
Hk  =  l. 


In  general,  this  approach  is  valid  for  three  dimensions.  However,  in  the  current 
context,  the  following  is  presented  for  the  two-dimensional  case  where  £l  «  {  and 
f3  =  T).  The  components  of  the  metric  tensor  may  then  be  written  as 

0u  =  *<  +  y?  (8?) 

013  =  x(x„  +  y(i b,  (88) 

022  =  +  (89) 

The  metric  tensor  can  be  decomposed  into  its  determinant  and  into  terms 

which  quantify  the  aspect  ratio  and  orthogonality  of  the  grid  cells.  The  skewness, 
or  the  angular  deviation  between  the  physical  and  computational  coordinate  lines, 
is  also  defined  via  components  of  the  metric  tensor.  In  terms  of  the  metric  tensor, 
the  inverse  Jacobian  is  defined  by 
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(j-1)3  =  detfa]  =  (x(y„  -  *„y<)3. 


The  orthogonality  of  the  grid  is  characterized  by  the  angle  9tJ  between  the  tangent 
lines  to  the  (  and  rj  coordinate  lines  at  the  points  of  intersection,  where 


-  Jk 

(i  ft  jr;  no  sum  on  ij) 

(90) 

+  y«SA, 

(911 

The  aspect  ratio  is  defined  by 

A>  =  M 

(k=2,n;  no  sum  on  k) 

(92) 

V 

_  + 

'J  *?+*?’ 

(93) 

while  the  cosine  of  the  skewness  angle,  a 

„  is  determined  by 

cos(a,) 


(94) 

(95) 


A  complete  error  analysis  of  the  Beam  and  Warming  algorithm  incorporating 
the  effects  of  the  metrics  is  beyond  the  scope  of  this  report  (see  McRae  and  Klopfer 
(12)).  However,  an  examination  of  the  metric  parameters  listed  above  can  provide 
some  insight  into  the  degree  of  truncation  error  induced  by  the  grid.  In  particular, 
one  should  suspect  increased  truncation  error  when  the  gradient  of  the  metrics  and 
the  associated  metric  parameters  becomes  large.  This  is  particularly  true  if  large 
flow  gradients  develop  at  the  sarnie  location. 
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III.  Numerical  Solution  of  the  NS  Equations 


The  numerical  procedure  for  the  solution  of  the  Navier-Stokes  equations  is  out¬ 
lined  in  this  chapter.  The  finite-difference  scheme  employed  is  based  on  the  work  by 
R.  M.  Beam  and  R.  F. Warming  and  is  known  as  the  Beam  and  Warming  approx¬ 
imate  factorization  algorithm  (5).  The  algorithm  has  four  pertinent  characteristics 
(19:139): 


1.  unconditional  stability, 

2.  spatial  factorization, 

3.  delta  formulation,  and 

4.  three-level  implicit  structure. 


S.l 


Implicit  NS  Code 


The  implicit  implementation  solves  the  strong  conservation  form  of  the  Navier- 
Stokes  equations.  In  delta  formulation,  with  first-order  Euler  time  differencing,  the 
scheme  may  be  written  as  follows  (18): 


d? 

* 


f  (dAn  d’Af’Nl  f  (dBn  trrn 

l  +  Uf  a?  +  a**)} 


-Al  |  -  Vl  -  V,f  +  -  W, 


drj  dr j3 

6_ 

<V 


~w,r\ 


(96) 


The  overall  scheme  progresses  in  time  as 


Un+ 1  =  Un  +  A17, 


(97) 


where  n  represents  the  temporal  index  and  A,  B,  M,  and  N  represent  the  Jacobian 
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matrices  (see  Appendix  B) 


A 

M 


dEi 

dU 

d}\ 

dU( 


(98) 

(99) 


Application  of  second-order  central  differencing  for  the  spatial  derivatives  yields  as 
a  function  of  nodal  position  (18:19) 


(100) 

-At  -  V2)tJ  +  nn(Ei  -  Wi),^} 

{/  +  A(  (^BtJ  -  tfKj)  }  A  0„  =  A  if*, 

(101) 

b = (« -  i)Ae 

1  <  i  <  IL 

(102) 

*7j  =  (>  “  l)Aq 

1  <  3  <  JL  , 

(103) 

with  the  finite-difference  operators  defined  as  (18:20) 


hUd  = 
Ww  = 

Mc/.j  « 

/*»>/»  J  ~ 


Af 

(J>j+  1/3  ~  ftj-i/ a) 

Aq 

(/»+ ij  ~  ft-ij) 
2A( 

(ftj+l  ~  /y-l) 
2AtJ 


(104) 

(105) 

(106) 
(107) 


The  transformation  derivatives  (metrics)  are  computed  directly  from  the  grid. 
Second-order,  centra]  differencing  is  employed  for  interior  points,  while  second-order, 
one-sided  differences  are  employed  on  the  boundaries.  The  scheme  is  implemented  in 


an  alternating-direction  fashion  by  first  solving  equation  (100)  along  lines  of  constant 
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(  for  2  <  j  <  JL-1  followed  by  the  solution  of  eqn  (101)  along  lines  of  constant  q  for 
2  <  *  <  IL-1.  This  results  in  a  block  tridiagonal  system  dong  each  coordinate  line 
due  to  the  computational  stencil  brought  about  by  the  approximate  factorization. 
The  boundary  equations  are  then  solved  explicitly. 

Compressible  flow  computations  are  suceptible  to  numerical  instabilities  in¬ 
duced  by  nonlinear  effects  due  to  large  flow  gradients  (e.g.  shock  waves  and  flow 
separation).  These  nonlinearities  are  especially  prone  to  producing  numerical  in¬ 
stabilities  at  high  Reynolds  numbers,  which  in  turn  can  slow  down  convergence 
considerably  (13:322). 

The  oscillations  induced  by  numerical  instabilities  can  be  damped  by  the  in¬ 
clusion  of  second-  and  fourth-order  artificial  viscosity  in  the  algorithm.  An  explicit 
fourth-order  damping  term,  Dt)  is  included  in  the  right-  hand  side  of  equation  (100). 
Two  second  order  damping  terms,  and  Dn  are  included  within  the  respective 
implicit  operators  in  equation  (101). 

The  equations  for  these  added  terms  are  (18:20) 


Dt  = 

(108) 

D(  = 

-u,A(  J.j/ 

(109) 

£ 

II 

-wAUjip'jI, 

(110) 

where  u>e  ~  1  and  w,  >  2wt. 

In  order  to  accelerate  convergence  for  a  steady  flow  computation,  a  local  time 
step,  Atli7,  is  introduced.  The  local  time  step  is  a  function  of  the  local  flow  con¬ 
dition  and  the  transformation  metrics.  The  time  step  is  controlled  by  the  user  via 
specifation  of  the  Courant  number,  C,  which  is  utilized  in  the  following  fashion 


The  maximum  local  timestep,  A as  mentioned,  is  then  a  function  of  the  local 
transformation  metrics  and  flow  conditions 


At 


A  S( 
A  Sn 


B  + 

M+av 

(A Sr3  +  A S’-*  +  M  +  e  ) 

Atj  ’ 

'  *  n  pAS*\Pr  Prt) 

■Jt\  +  S^A( 


(112) 

(113) 

(114) 


Generally,  for  an  explicit  algorithm,  the  Courant  number  must  be  less  than 
one  for  stability.  However,  for  the  Beam  and  Wanning  implicit  algorithm  it  can  be 
much  higher  and  is  typically  set  to  twenty  (5:20). 


3.2  Convergence  Criteria 

By  considering  only  steady-state  solutions,  the  specification  of  convergence 
criteria  is  simplify  d.  In  general,  some  small  oscillation  in  flow  variables  will  persist 
indefinitely.  These  oscillations  may  be  too  small  to  appreciably  affect  the  solution. 
A  suitable  convergence  criteria  must  be  established  to  ascertain  when  the  oscillations 
axe  satisfactorily  small. 

The  convergence  criteria  used  in  this  study  follows  the  example  of  Visbal  (18) 
and  Boyles  (5).  Convergence  was  assesed  by  monitoring  the  airfoil  lift  and  drag  coef¬ 
ficients.  The  solution  was  assumed  converged  when  the  amplitude  of  the  oscillations 
for  Ct  was  less  than  0.1  percent  and  the  variation  in  Cj  was  less  than  one  drag  count, 
where  one  drag  count  equals  0.0001. 


IV.  Results  and  Discussion 


This  chapter  presents  the  results  obtained  from  the  numerical  study  of  the  two 
differing  grid  topologies.  The  study  is  divided  into  two  basic  parts.  The  first  consists 
of  a  grid  refinement  study  using  the  C-grid.  Progressively  finer  grids  are  employed 
for  subsonic  and  tram  sonic  cases.  The  second  part  is  the  comparative  study  between 
the  C-  and  O-grid  topologies.  An  overview  of  the  various  grids  employed  and  their 
construction  is  presented  first.  This  is  followed  by  the  results  of  the  C-  grid  refinement 
study  and  a  comparative  analysis  between  the  C-  and  O-grid. 


1  Grids  Employed 

There  are  procedural  differences  in  the  construction  of  the  two  grid  types.  The 
construction  of  C-grids  is  based  on  a  straightforward  application  of  the  Barth  and 
Kinsey  code  (9).  The  construction  of  O-grids  is  complicated  by  the  need  to  enforce  a 
finite  curvature  at  the  trailing  edge,  allowing  the  O-grid  generator  to  function.  The 
first  step  is  the  construction  of  the  airfoil  shape.  The  airfoil  is  then  truncated  at 
a  specified  point  and  a  circular  arc  is  affixed  to  form  the  trailing  edge.  The  point 
distribution  around  the  airfoil  is  then  specified  and  the  actual  generation  of  the  grid 
can  then  proceed.  This  procedure  represents  an  added  degree  of  difficulty  which 
should  be  considered  when  comparing  the  two  types  of  grid. 

The  grids  employed  in  this  study  are  listed  in  Tables  1  and  2.  Figures  2 
through  5  display  the  progressively  finer  C-grids  employed.  C-Grids  5  through  8, 
shown  in  Figures  6  through  9  were  constructed  for  specialized  purposes.  C-grid 
5  is  an  asymmetric  grid,  constructed  with  more  points  and  finer  spacing  on  the 
upper  surface  to  test  for  better  shock  capturing  in  the  transonic  case.  C-grid  6  was 
constructed  with  a  finer  leading  edge  spacing  in  am  attempt  to  better  match  the 
experimental  pressure  coefficient  data  at  that  location.  C-grid  7  was  constructed  in 
an  attempt  to  smooth  the  gradient  in  the  tramsformation  metrics,  especially  at  the 
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trailing  edge/ wake  junction.  C-Grid  8  represents  a  further  attempt  to  smooth  the 
metrics  without  sacrificing  fine  tangential  spacing  at  the  trailing  edge.  The  three 
basic  O-grids  shown  in  Figures  10  through  12  vary  at  the  point  in  which  the  nominal 
airfoil  was  truncated,  thus  they  vary  in  the  radius  of  curvature  at  the  respective 
trailing  edges.  O-grid  4,  shown  in  Figure  13,  was  constructed  with  a  finer  grid 
spacing  on  the  upper  surface  to  test  for  better  shock  resolution.  Figure  14  shows  an 
example  of  the  trailing-edge  geometry  for  the  O-grids,  while  Figures  15  and  16  offer 
a  comparison  of  leading-edge  geometries  for  the  C-  and  O-grids. 

4- 2  C-Grid  Refinement  Study 

Two  specific  flow  regimes  were  considered  for  the  C-grid  refinement  study;  a 
subsonic  and  transonic  case.  The  subsonic  case  chosen  was  at  M  =  0.3,  a  =  4.04 
degrees,  and  Re  =  1,860,000.  The  transonic  case  chosen  was  at  M  =  0.775,  a  =  2.03 
degrees,  and  Re  =  1,900,000.  These  cases  were  chosen  because  of  the  availability  of 
comparative  experimental  data  (2)  (11). 

Figures  17  and  18  display  the  variation  of  pressure  coefficient  on  the  surface 
for  the  subsonic  case.  The  differences  in  the  plots  appear  most  readily  at  the  leading 
edge.  C-grids  1  and  4  display  almost  identical  results,  each  displaying  a  spike  in  Cp  on 
the  upper  surface  near  the  leading  edge  and  generally  overestimating  the  magnitude 
of  Cp  along  the  upper  surface  to  the  midchord  point.  A  close  examination  reveals  that 
C-grid  1  matches  the  experimental  data  slightly  better  than  C-grid  4  in  that  range. 
However,  an  examination  of  similar  experimental  data  reveals  the  sensitivity  of  the 
pressure  distribution  to  angle  of  attack  (7:58).  The  corrected,  experimental  angle 
of  attack  (as  opposed  to  the  geometric  angle)  is  a  difficult  parameter  to  determine 
precisely,  and  is  a  likely  cause  of  this  discrepancy.  The  results  using  a  finer  leading- 
edge  spacing  (C-grid  6)  are  shown  in  Figure  19.  The  spike  in  pressure  at  the  leading 
edge  is  slightly  reduced,  but  the  trend  remains  the  same.  The  numerical  results 
obtained  are  duplicated  using  a  finer  grid  (299x100)  by  Boyles  (5:43). 
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Table  1.  C-Grid  Specifications 


O-Grid 

... 

ILx  JL 

K2ZBI 

Pt»lo 

TEC-1 

Aile 

mm 

KDIST 

1 

204x108 

102 

102 

.0051c 

.0050 

.000150 

.0004 

15 

2 

204x108 
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102 

.0025c 

.0050 

.000100 

.0004 

15 

3 

204x108 
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.0014c 

.0050 

.000075 

.0004 

15 

4 

204x108 

102* 
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.0014c 
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.000075 
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15 

*  Grid  spacing  refined  along  forward  surface 


Table  2.  O-Grid  Specifications 
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The  results  for  the  transonic  case  show  a  greater  variation  with  grid  spacing. 
Figures  20  and  21  display  a  convergence  toward  the  experimental  data  with  grid 
refinement.  As  the  grid  becomes  finer  and  the  chordwise  distribution  of  points  in¬ 
creases,  the  numerical  solution  approaches  the  experimental  data.  C-Grid  5  was 
constructed  in  an  attempt  to  fit  the  grid  to  the  expected  solution.  The  number  of 
points  on  the  upper  surface  was  increased  and  the  chordwise  spacing  decreased  at 
the  expected  location  of  the  shock.  A  slight  variation  of  the  pressure  distribution  is 
evident  in  Figure  22,  with  the  shock  definition  slightly  enhanced. 

In  general,  there  will  be  some  variation  from  experimental  data  due  to  a  com¬ 
bination  of  experimental  error,  numerical  truncation  error,  and  the  need  to  apply 
artificial  viscosity  to  achieve  solution  convergence.  To  determine  the  effect  of  the 
second-order  artificial  viscosity,  the  value  utilized  in  the  previous  case  was  reduced 
by  fifty  percent.  Figure  23  displays  the  result  using  C-grid  4;  the  reduction  in 
artificial  viscosity  produced  an  even  closer  agreement  with  the  experimental  data. 
Comparing  Figures  22  and  23  indicates  that  minimizing  the  artificial  viscosity  may 
be  as  or  even  more  important  in  obtaining  an  accurate  solution  than  am  overreliauice 
on  reducing  grid  spacing.  The  importance  of  minimizing  artificial  viscosity  is  further 
demonstrated  by  Figure  24,  in  which  the  effect  of  axtificial  viscosity  on  Cj  and  C*  is 
shown.  In  particular,  the  value  of  Ci  varies  by  nearly  10  percent,  while  C\  varies  by 
nearly  4  percent  over  the  range  of  artifical  viscosities  considered. 

The  contour  plots  display  the  same  trends  as  the  plots  of  Cp  on  the  surface. 
Comparison  of  the  Mach  contours  in  Figures  25  and  26  reveals  that  the  solution  ob¬ 
tained  with  C-grid  1  exhibits  a  misplaced  sonic  line,  moves  the  shock  too  fax  forward, 
and  has  a  poor  resolution  of  the  shock  wave.  Figure  27  displays  the  experimental 
location  of  the  shock  at  the  midchord  point  (11:69).  The  solution  obtained  using 
C-grid  4  places  the  shock  at  the  correct  location  amd  displays  better  shock  resolution. 

Figures  28  and  29  exhibit  the  dependency  of  convergence  rate  on  flow  conditions 
and  grid  characteristics.  The  plots  display  the  oscillation  of  lift  coefficient  with 


iteration  number,  N.  Convergence  is  obtained  (or  the  subsonic  case  (M  =  0.3)  in  less 
than  600  iterations  using  C-grid  2.  However,  increasing  the  Mach  number  to  0.775 
produces  a  marked  oscillation  in  Ct  which  persists  beyond  5000  iterations.  This  is 
attributed  at  least  in  part  to  the  point  distribution  in  the  normal  direction.  Referring 
to  Table  1,  C-grid  2  uses  only  50  percent  more  points  in  the  normal  direction  to  cover 
twice  the  distance  (i.e.,  ten  chord  lengths  instead  of  five).  The  poor  point  distribution 
can  contribute  to  numerical  error,  which  is  then  further  exacerbated  by  large  flow 
gradients.  When  additional  points  are  added  in  the  normal  direction,  the  problem 
ceases,  as  revealed  by  Figure  30. 


4-S  Comparison  of  C-grid  and  O-Grid  Results 

Three  regimes  are  chosen  for  comparative  purposes:  subsonic,  transonic,  and 
supersonic.  Experimental  data  is  used  as  the  yardstick  of  measure  in  the  subsonic  and 
transonic  cases,  while  the  differences  in  flow  structure  is  considered  in  the  supersonic 
case.  In  general,  the  higher  the  Mach  number,  the  more  differences  appear  in  the 
solutions  between  the  two  grid  types.  Only  small  differences  are  noticed  in  the 
subsonic  regime,  while  more  significant  differences  are  noted  for  the  other  cases. 


4.S.J  Subsonic  Case  Of  the  first  three  O-grids  employed,  (see  Table  2)  the 
only  difference  is  the  point  at  which  the  trailing  edge  is  truncated  and  the  circular 
arc  affixed.  Of  course,  the  further  aft  this  truncation  occurs,  the  smaller  the  radius 
of  curvature  of  the  arc.  In  general  there  was  little  difference  in  the  solutions  obtained 
using  these  three  grids.  Comparing  Figures  31  and  32  reveals  no  discernable  differ¬ 
ence  in  the  plots.  Figure  31  reveals  the  effect  of  the  truncation  at  approximately 

1  =  .97,  where  the  data  is  similarly  truncated.  However,  the  pressure  distribution 
over  the  remainder  of  the  surface  is  largely  unaffected.  The  integrated  variables, 
CiandCi,  are  also  comparable.  The  variation  of  C|  and  C*  is  less  them  2  percent 
when  comparing  O-grids  1  and  3,  and  less  than  0.5  percent  when  comparing  O-grids 

2  and  3. 
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Comparing  these  results  with  those  of  C-grids  4  and  6  for  the  subsonic  case 
(Figures  18  and  19)  shows  a  similar  agreement  with  two  exceptions:  the  leading 
and  trailing  edges  show  slight  differences.  At  the  leading  edge,  C-grid  4  displays  a 
greater  pressure  spike  on  the  upper  surface.  Only  when  the  leading-edge  spacing  is 
reduced  to  .003  (C-grid  5)  do  the  pressure  values  agree  very  closely  at  that  location. 
When  equal  spacing  is  specified  at  the  leading  edge,  the  O-grid  performs  slightly 
better.  Note  that  the  specification  of  equal  spacing  does  not  enforce  an  equal  point 
distribution  in  the  region  near  the  leading  edge  for  the  C-  and  O-grids.  This  can  be 
seen  in  Figures  5  and  12;  the  point  distributions  are  similar  but  not  exact.  There 
is  also  some  variation  in  Cp  at  the  trailing  edge.  The  averaging  procedure  applied 
to  the  boundary  conditions  in  the  wake  for  the  case  of  the  C-Grid  may  induce 
some  oscillatory  behavior  there.  The  O-grid  has  no  branch  cut  in  the  wake  and  so 
this  difficulty  does  not  arise.  See  Appendix  C  for  a  discussion  of  these  boundary 
conditions. 


The  comparison  of  lift  coefficient  versus  angle  of  attack  is  presented  in  Fig¬ 
ure  33.  Results  are  presented  for  Reynolds  number  6,000,000.  The  experimental 
data  (1:462)  represents  an  average  of  the  results  presented  for  this  value  of  Reynolds 
number.  The  agreement  with  this  data  indicates  that  the  turbulence  level  for  the 
experiment  could  not  be  precisely  matched,  since  the  only  differences  in  the  experi¬ 
mental  data  was  the  transition  point  and  turbulence  level.  Transition  was  specified 
at  the  leading  edge  for  both  C-  and  O-grids  in  the  numerical  code,  attempting  to 
match  the  experimental  case  where  .011  inch  gTain  was  applied  at  the  leading  edge. 
The  numerical  results  show  excellent  agreement  up  to  10  degrees  angle  of  attack. 
Past  this  point  there  is  some  variation,  with  the  O-grid  correctly  predicting  artaii  at 
approximately  12  degrees.  The  data  curve  truncates  at  the  twelve  degree  point  for 
the  C-grid  because  the  solution  diverged  at  higher  values. 


4-S.2  Transonic  Case  In  the  transonic  case,  there  again  is  no  discernable 
difference  in  the  pressure  distribution  on  the  surface  between  O-grids  1,  2,  and  3, 
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as  shown  by  Figures  34  and  35.  The  data  is  truncated  near  the  trailing  edge  for  0- 
grid  3,  as  noted  in  the  subsonic  case,  but  the  overall  distribution  remains  unaffected. 
However,  a  comparison  with  the  result  of  C-grid  4,  also  shown  Figure  34,  indicates 
that  the  C-grid,  with  equal  values  set  for  second-  and  fourth-order  artificial  viscosity, 
performed  slightly  better  in  predicting  Cp  in  the  vicinity  of  the  shock  wave.  In  an 
attempt  to  match  this  result,  the  chordwise  grid  spacing  was  reduced  forward  of 
the  shock  wave  for  O-grid  5.  However,  the  C-grid  still  outperformed  the  O-grid 
in  matching  the  experimental  data,  as  the  results  of  Figure  36  fail  to  show  any 
improvement. 

A  comparison  of  Mach  contours  shows  two  deficiencies  of  O-grid  3  vis-4-vis 
C-grid  4.  See  Figures  25  and  37.  Firstly,  the  O-grid  displays  a  poorer  resolution 
of  the  shock  wave,  with  the  shock  wave  positioned  too  far  forward.  Secondly,  the 
O-grid  displays  large  errors  in  the  wake  region.  The  O-grid  has  a  large  concentration 
of  grid  points  directly  aft  of  the  trailing  edge.  However,  beyond  x  w  1.3  the  grid 
spacing  becomes  larger.  This  inadequate  resolution  produces  the  choppy,  stairstep 
type  of  contour  seen  in  the  wake  of  Figure  37.  O-grid  5  was  constructed  specifically  to 
address  the  first  problem.  The  result  is  seen  in  Figure  38.  The  resolution  of  the  shock 
wave  for  this  case  equals  that  produced  by  the  C-grid.  In  addition  the  small  contour 
directly  aft  of  the  shock  wave  is  nearly  identical  to  that  of  the  C-grid,  in  contrast  to 
the  previous  case.  The  O-grid  solution  still  places  the  shock  wave  slightly  forward 
of  the  C-grid  shock  wave,  but  the  overall  definition  is  much  improved.  However,  the 
shifting  of  points  to  the  surface  has  the  result  of  reducing  still  further  the  number 
of  points  in  the  wake.  The  deficiency  noted  in  the  wake  for  O-grid  3  thus  becomes 
even  more  severe  for  O-grid  4. 

Interestingly,  contours  of  Cp  indicate  a  deficiency  of  the  C-grid.  Upon  exam- 
intion  of  Figure  39,  the  contour  in  the  wake  reveals  a  small  oscillation  that  occurs 
directly  along  the  branch  cut.  This  oscillation  is  most  likely  due  to  the  handling  of 
the  boundary  conditions  along  the  cut.  The  corresponding  O-grid  plot,  Figure  40, 


has  a  similar  contour  in  the  wake,  but  lacks  the  oscillation  apparent  in  the  C-grid 
case.  In  general,  the  same  comments  apply  to  the  pressure  contours  in  terms  of  shock 
location  and  resolution  as  for  the  Mach  contours,  i.e.,  the  O-grid  solutions  display 
comparable  shock  resolution  only  if  a  particular  effort  is  made  to  concentrate  points 
in  the  vicinity  of  the  shock. 

The  comparison  of  vorticity  between  the  two  grids,  shown  in  Figures  42  and 
43,  display  similar  profiles.  However,  the  inadequate  resolution  of  the  O-grid  wake 
is  again  apparent  in  the  step-like  curves  occuring  there. 

4.S.3  Supersonic  Case  The  supersonic  case  considered  is  not  compared  with 
experimental  data  as  in  the  previous  two  cases.  The  conditions,  M  =  1.1,  a  =  6°, 
exceed  the  design  specifications  of  the  airfoil  and  experimental  data  in  this  regime  is 
very  limited  and  somewhat  suspect  (10:3).  However,  a  comparison  can  still  be  made 
based  on  global  flow  characteristics. 

Contour  plots  of  Mach  contours,  Cp,  and  vorticity  are  presented  in  Figures  44 
through  51.  All  of  the  O-grid  plots  display  the  weakness  noted  earlier  of  inadequate 
resolution  in  the  wake.  The  contours  there  are  generally  choppy  and  ill-formed,  just 
as  in  the  previous  cases.  The  resolution  of  the  lambda  shock  forming  at  the  trailing 
edge  is  generally  quite  good  for  the  O-grid,  superior  to  that  of  the  C-grid.  This  is 
apparent  for  both  the  Mach  and  Cp  contours.  An  examination  of  O-grid  3  (Figure 
12)  reveals  why  this  is  the  case;  the  O-grid  typically  has  a  pronounced  clustering 
of  points  near  the  trailing  edge,  nearly  coincident  with  the  formation  of  the  shock 
wave.  This  aids  the  resolution  of  the  shock  wave,  at  least  near  the  body. 

The  corresponding  C-grid  plots  (Figures  45,  47,  and  49)  display  an  obvious 
shortcoming  in  the  large  amounts  of  numerical  noise  in  the  region  around  the  shock 
wave.  These  disturbances  ate  certainly  nonphysical  and  appear  most  strongly  just 
prior  to  the  shock  wave.  They  appear  to  dampen  as  they  propogate  through  the 
shock  and  further  downstream.  The  plot  of  vorticity  is  most  interesting,  as  the 


disturbances  appear  to  outline  the  grid  elements  themselves.  It  can  be  concluded 
that  these  disturbances  do  not  arise  from  the  handling  of  the  branch  cut  in  the 
wake,  previously  implicated  in  other  cases  for  two  reasons.  Firstly,  the  flow  remains 
supersonic  after  the  oblique  shock  and  so  downstream  errors  should  not  propogate 
upstream,  and  secondly,  the  disturbances  appear  to  dampen  out,  not  magnifying,  as 
they  approach  the  wake.  For  these  reasons,  the  most  likely  source  of  this  error  is 
the  transformation  metrics  themselves.  This  hypothesis  was  tested  by  construcing 
C-grid  8  in  an  attempt  to  smooth  the  metrics  in  the  region  where  the  numerical 
noise  was  generated.  The  results  obtained  utilizing  C-grid  8  are  shown  in  Figures  50 
and  51.  The  Mach  contours  display  very  little  numerical  noise  in  the  region  forward 
of  the  shock,  however,  numerical  noise  persists  aft  of  the  shock  wave.  The  vorticity 
contours  display  the  same  trend.  The  large  amounts  of  numerical  noise  forward  of 
the  shock  wave  evident  in  Figure  49  are  completely  absent  in  Figure  51.  However, 
the  numerical  noise  persists  aft  of  the  shock  wave.  See  the  section  on  the  comparison 
of  the  transformation  metrics  for  a  further  discussion  of  this  point. 

4.8.4  Convergence  Rate  Comparrison  There  are  minor  differences  in  the  con¬ 
vergence  properties  for  the  C-  and  O-grids.  Plots  of  Cj  and  C*  are  presented  in  Fig¬ 
ures  52  through  55  for  the  subsonic  and  transonic  cases.  In  general,  the  oscillations 
in  Ci  and  Cj  are  more  pronounced  for  the  O-grid  than  the  C-Grid  for  the  subsonic 
case.  Convergence  is  attained  for  the  O-Grid  solution  in  about  1400  iterations,  while 
the  oscillations  occuring  in  the  case  of  the  C-grid  exhibit  less  amplitude,  with  conver¬ 
gence  attained  in  about  1100  iterations.  For  the  transonic  case,  the  C-grid  solution 
converges  in  1700  iterations  while  the  O-grid  solution  converges  in  2000  iterations. 

This  difference  in  convergence  rate  is  due  in  part  to  the  smaller  minimum  time 
step  which  must  be  specified  to  avoid  divergence  in  the  case  of  the  O-grid.  The  local 
time  step,  being  a  function  of  grid  spacing,  essentially  becomes  too  large  around 
the  trailing  edge  of  the  O-grid,  where  the  tangential  spacing  is  extremely  small. 
Decreasing  the  user  specified  time  step  with  respect  to  the  C-grid  value  allows  the 


Q-grid  solution  to  converge.  It  is  significant  that  there  are  particular  cases  for  which 
the  O-grid  solution  converges  while  the  C-grid  solution  diverges.  For  example  at 
M  —  0.3,  a  —  13°,  and  Re  =  6,000,000  the  O-grid  solution  converges  while  the 
C-grid  solution  diverges.  This  occurs  despite  lowering  the  minimum  timestep  for 
the  C-grid  (as  well  as  the  CFL  input).  The  typical  mechanism  for  divergence  is  the 
appearance  of  a  negative  value  of  temperature  in  the  field.  In  terms  of  convergence 
criterea,  the  C-grid  solution  exhibits  less  oscillation  in  the  integrated  variables  and 
tends  to  converge  somewhat  sooner,  while  the  O-grid  solution  converges  for  more 
critical  flow  conditions  when  the  C-Grid  solution  diverges.  This  suggests  that  a 
different  source  of  numerical  error  is  occuring  to  produce  the  difficulties  mentioned. 

4-4  Comparison  of  Transformation  Metrics 

There  is  a  significant  difference  in  the  appearance  of  the  transformation  metrics 
between  the  C-grid  and  O-grid  topologies.  It  is  most  intuitive  to  display  the  metrics 
in  terms  of  geometric  parameters  such  as  the  Jacobian,  orthogonality,  and  aspect 
ratio  (8).  Figures  56  and  57  display  aspect  ratio  contours  for  C-grid  4  and  O-grid  3 
respectively.  The  basic  structure  of  the  contours  is  similar,  there  are  some  differences 
in  the  leading-edge  area,  while  both  display  a  steeper  gradient  in  the  area  above  and 
below  the  trailing  edge.  The  contours  for  the  O-grid  are  generally  more  rounded, 
while  the  C-grid  displays  a  cusp  in  the  transition  to  the  wake  region. 

The  departure  from  orthogonality  is  displayed  in  Figures  58  and  59.  Both 
grids  are  nearly  orthogonal  in  the  wake,  i.e.,  coj(0)  «  0.  The  largest  differences  here 
are  at  the  leading  edge.  Surprisingly,  the  C-grid  exhibits  a  larger  departure  from 
orthogonality  in  the  region  around  the  leading  edge  and  above  and  below  the  leading 
edge.  This  may  partially  explain  why  the  C-grid  requires  a  finer  grid  spacing  at  the 
leading  edge  to  match  the  more  accurate  O-grid  performance  there,  as  discussed 
earlier. 

Contour  plots  of  the  Jacobian,  or  the  area  ratio  between  the  physical  and 


computational  domains,  are  presented  in  Figures  60  through  63.  Some  significant 
differences  are  evident  between  the  C-  and  O-grid  cases.  The  C-grid  displays  a 
distinct  cusp  above  and  below  the  trailing  edge  of  the  airfoil.  The  O-grid  has  a  much 
smoother  transition  in  this  area  and  the  contours  display  such  cusp.  The  numerical 
error  will  be  a  function,  in  part,  of  the  gradient  in  the  transformation  metrics. 
Therefore,  at  the  cusp  point,  the  C-grid  will  generate  metric-induced  errors  that  are 
larger  than  those  of  the  O-grid.  This  offers  a  possible  explanation  for  the  numerical 
oscillations  occuring  in  the  supersonic  case.  The  oscillations  occur  in  precisely  the 
region  where  the  gradient  of  the  Jacobian  displays  a  discontinuity,  corresponding  to 
a  cusp  in  the  Jacobian  contours.  The  O-grid  displays  a  well  behaved  gradient  there, 
and  no  such  oscillations  occur. 

In  an  attempt  to  rectify  this  situation,  C-grid  7  was  generated  with  the  intent 
of  smoothing  the  metric  gradients  (see  Figure  8).  It  was  hoped  that  this  would  al¬ 
low  a  better  solution  and  perhaps  allow  convergence  for  several  C-grid  cases  which 
diverged.  However,  this  attempt  was  a  failure  because  the  means  chosen  to  accom¬ 
plish  the  smoothing  was  inadequate.  Smoothing  was  accomplished  by  transfering 
points  from  the  airfoil  surface  to  the  wake  and  increasing  the  specified  spacing  at 
the  trailing  edge.  Figure  62  reveals  that  significant  smoothing  of  the  Jacobian  con¬ 
tours  was  achieved.  However,  this  apparently  increased  the  total  truncation  error 
due  to  increased  grid  spacing  to  the  extent  that  any  improvement  due  to  smoothing 
the  metrics  was  not  discernable.  A  difficult  trade-off  occurs  within  the  constraint 
of  a  fixed  number  of  grid  points.  This  involves  allowing  sufficient  grid  resolution  in 
critical  flow  areas  while  dispermiting  any  abrupt  changes  in  the  metrics.  Optimizing 
this  trade  off  may  reduce  the  total  truncation  error. 

C-grid  8  (see  Figure  9)  was  also  generated  with  the  intent  of  smoothing  the 
metric  gradients.  The  goal  was  to  reduce  the  large  amounts  of  numerical  noise 
evident  in  the  C-grid  solution  for  the  supersonic  case  (see  Figures  45  and  49).  The 
method  employed  to  achieve  smoothing  was  to  enforce  uniform  tangential  spacing 


from  x  »  0.8  to  x  =  1.0  (the  trailing  edge).  The  result  on  the  Jacobian  contours  is 
evident  in  Figure  63.  The  smoothing  is  effective  forward  of  the  trailing  edge,  but  the 
cusp  still  remains  at  x  =  1.0,  since  the  uniform  spacing  was  not  maintained  beyond 
this  point.  There  is  a  strong  correlation  between  the  Jacobian  contours  and  the 
vorticity  contours  seen  in  Figure  51.  The  large  amounts  of  numerical  noise  evident 
in  Figure  49  is  completely  gone  in  the  region  forward  of  the  shock  wave  in  Figure 
51.  This  corresponds  to  the  region  in  Figure  63  where  the  metric  smoothing  occurs. 
However,  the  numerical  noise  persists  aft  of  the  shock  wave,  since  the  magnitude 
of  the  metric  gradients  are  unaltered  beyond  x  -  1.0.  The  overall  reduction  in  the 
level  of  numerical  noise  could  also  be  attributed  to  finer  grid  spacing  (irregardless 
of  the  metrics).  However,  a  comparison  with  the  midchord  shock  wave  occuring 
in  the  transonic  case  (see  Figure  26),  shows  no  numerical  noise  evident,  despite  a 
relatively  course  tangential  spacing.  It  is  pertinent  that  the  Jacobian  contours  are 
well  behaved  in  the  midchord  region.  Therefore,  it  is  inferred  that  the  coincidence  of 
large  flow  gradients  and  large  metric  gradients  contribute  to  the  type  of  numerical 
noise  seen  in  Figure  49. 


Figure  2.  C-grid  1  (119x40) 
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Figure  5.  C-grid  4  (219x100) 
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Figure  6.  C-gnd  5:  Refined  Spacing  for  Shock  (229x100) 
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Figure  7.  C-grid  6:  Modified  Leading-Edge  Spacing  (219x100) 


Figure  8.  C-grid  7:  Metrics  Smoothed  (219x100) 
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Figure  9.  C-grid  8:  Metrics  Smoothed  (299x60) 
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Figure  10.  O-grid  1:  .0051c  TEC  (204x108) 
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Figure  11.  O-grid  2:  .0025c  TEC  (204x108) 
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Figure  12.  O-grid  3:  .0015c  TEC  (204x100) 
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Figure  13.  O-grid  4:  ,0015c  TEC  Modified  Spacing  (204x108) 
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Figure  14.  O-grid  2:  Trailing-Edge  Geometry 
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Figure  15.  O-grid  3:  Leading  Bdge  Geometry 
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Figure  17.  C-grid  1:  Subsonic  Pressure  Coefficient  Profile 


50 


Figure  18.  C-grid  4:  Subsonic  Pressure  Coefficient  Profile 
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Figure  19.  C-grid  5:  Subsonic  Pressure  Coefficient  Profile 


Figure  20.  C-grid  1:  Transonic  Pressure  Coefficient  Profile 
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Figure  21.  C-grid  4:  Transonic  Pressure  Coefficient  Profile 
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Figure  22.  C-grid  5:  Transonic  Pressure  Coefficient  Profile 


Figure  23.  Transonic  Pressure  Coefficient  Profile,  Reduced  Damping 
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Second-Order  Artificial  Viscosity 

Figure  24.  Ct  and  C*  ?s  Artificial  Viscosity 
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Figure  27.  Experimental  Location  of  Shock  Wave,  M  •  0.775,  a  -  2.03 

(Ref:  U) 
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Figure  28.  C-grid  2:  Lift  Coefficient  vs  Iterations,  M  s=  0.3,  a  =  4.04 


6 


N 


Figure  30.  C-grid  4:  Lift  Coefficient  vs  Iterations,  M  =  0.775,  a  =  2.03 
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Figure  32.  O-grid  3:  Subsonic  Pressure  Coefficient  Profile 
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Figure  34.  C-grid  4,  O-grid  3:  Transonic  Pressure  Coefficient  Profiles 
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Figure  35.  O-grid  2:  Transonic  Pressure  Coefficient  Profile 
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Figure  36.  O-grid  4:  Transonic  Pressure  Coefficient  Profile 
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Figure  39.  C-grid  4:  Pressure  Coefficient  Contours,  M  =  0.775,  a  —  2.03 
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Figure  46.  O-grid  3:  Pressure  Coefficient  Contours,  M  =  1.1 
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Figure  60.  C-grid  4:  Jacobian  Contours 
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V.  Conclusions  and.  Recommendations 


The  numerical  solutions  obtained  with  both  the  C-  and  O-grid  configurations 
using  the  Beam  and  Warming  approximate  factorization  algorithm  compare  favor¬ 
ably  to  both  experiment  and  previous  numerical  results.  The  numerical  results  pro¬ 
duced  indicate  that  both  grids  are  satisfactory  in  a  general  application.  However, 
specific  weaknesses  arise  for  each  type  of  grid  in  certain  situations.  In  these  instances 
there  may  be  an  advantage  in  selecting  a  particular  type  of  grid. 

Studies  in  the  subsonic  range  indicate  that  when  equal  spacing  of  grid  points 
is  specified  at  the  leading  edge,  the  O-grid  is  superior  to  the  C-grid  in  determining 
the  pressure  coefficient  in  the  vicinity  of  the  leading  edge.  The  O-grid  also  displays  a 
pronounced  advantage  when  large  flow  gradients  appear  in  the  vicinity  of  the  trailing 
edge,  such  as  the  formation  of  a  lambda  shock.  The  O-grid  typically  exhibits  a  clus¬ 
tering  of  grid  points  around  the  trailing  edge,  and  the  gradient  of  the  transformation 
metrics  is  considerably  smoother.  For  certain  critical  flow  conditions,  e.g.,  a  combi¬ 
nation  of  high  Reynolds  number,  high  Mach  number  or  angle  of  attack,  the  O-grid 
displays  superior  numerical  stability  when  compared  with  the  C-grid.  This  result  is 
a  likely  result  of  two  effects:  the  global  smoothness  of  the  transformation  metrics 
compared  to  those  of  the  C-grid,  and  the  errors  induced  by  the  boundary  conditions 
across  the  C-grid  branch  cut.  The  increased  numerical  stability,  demonstrated  by 
the  O-grid,  allows  a  better  estimation  of  a  at  C\mmw. 

The  C-grid  displayed  several  advantages  over  the  O-grid  as  well.  One  advantage 
is  ease  of  application.  It  is  generally  easier  to  construct  a  C-grid  to  a  given  set 
of  specifications  than  an  O-grid.  The  C-grid  can  be  constructed  in  a  single  step 
with  minimal  effort,  while  the  O-grid  requires  multiple  steps.  It  is  also  difficult  in 
the  case  of  the  O-grid  to  concentrate  grid  points  in  the  wake  while  maintaining  a 
reasonable  spacing  at  the  trailing  edge.  The  O-grid  generator  tends  to  coalesce  points 
in  a  fan-like  shape  directly  above  and  below  the  trailing  edge.  Too  fine  a  spacing 
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at  the  trailing  edge  can  produce  numerical  instabilities  when  using  a  hyperbolic 
grid  generator.  Therefore,  it  becomes  difficult  to  maintain  adequate  grid  spacing 
in  the  wake.  In  the  wake  region,  both  grids  display  problems.  The  O-grid  has 
difficulty  resolving  flow  structures  because  of  the  large  grid  spacing  there.  The  C- 
grid  develops  small  oscillations  along  the  branch  cut  where  the  boundary  conditions 
are  applied.  However,  the  C-grid  handles  the  flow  resolution  better  in  the  wake  when 
the  flow  gradients  are  concentrated  directly  behind  the  airfoil.  In  this  region  the  C- 
grid  displays  a  considerably  higher  concentration  of  grid  points.  The  C-grid  also 
predicts  the  pressure  distribution  and  shock  location  more  accurately  the  O-grid  in 
the  transonic  case.  The  O-grid  can  approach  the  C-grid  performance  in  this  regard 
only  if  a  deliberate  attempt  is  made  to  match  the  point  distribution  occuring  in  the 
C-grid  case,  at  the  cost  of  a  further  deterioration  of  the  solution  in  the  wake. 

The  selection  of  a  specific  grid  type  should  be  geared  to  the  application  consid¬ 
ered,  as  indicated  by  the  conclusions  enumerated  above.  However,  there  are  several 
recommendations  that  could  increase  the  utility  of  both  types  of  grids.  In  the  case 
of  the  C-grid,  central  differencing  across  the  branch  cut  could  reduce  or  eliminate 
some  of  the  spurious  results  demonstrated  in  the  wake.  However,  this  would  create 
additional  computational  and  memory  requirements.  In  the  case  of  the  O-grid,  the 
lack  of  an  ability  to  control  grid  spacing  in  the  wake,  when  applying  a  hyperbolic 
grid  generator,  presents  a  significant  problem.  Therefore,  an  elliptic  grid  generation 
scheme  should  be  considered.  This  represents  an  added  degree  of  difficulty  due  to  the 
iterative  numerical  process  required  and  the  constraint  of  the  outer  computational 
domain.  However,  the  constraint  of  the  outer  domain  offers  a  means  to  control  the 
grid  spacing  in  the  wake. 

An  important  consideration,  which  applies  to  both  C-  and  O-grids,  is  the 
minimization  of  artificial  viscosity.  In  particular,  second-order  artificial  viscosity  has 
a  considerable  effect  upon  the  pressure  distribution  on  the  surface  and  the  integrated 
variables  Ct  and  C<t.  The  minimization  of  artificial  viscosity  should  be  considered  a 


prerequisite  to  obtaining  an  accurate  solution. 

Another  approach  worth  consideration  is  an  attempt  to  reduce  the  overall 
truncation  error  by  smoothing  the  transformation  metrics.  This  involves  a  trade  off 
between  the  need  to  concentrate  grid  points  in  an  area  of  interest  and  not  allowing 
a  steep  gradient  in  the  transformation  metrics  to  result.  The  most  likely  area  of 
application  would  be  in  the  region  above  and  below  the  trailing  edge  in  the  case  of 
the  C-grid.  A  method  of  smoothing  the  metrics  in  this  region,  without  unduly  com¬ 
promising  the  grid  spacing,  may  provide  a  means  of  reducing  the  overall  truncation 
error. 
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Appendix  A.  Non- dimensional  Variables 


The  Navier-Stokes  equations  axe  generally  nondimensionalized  to  allow  for  ease 
of  application.  This  is  accomplished  by  selecting  reference  length,  velocity,  and 
mass  scales,  L*,  {/£,,  and  respectively.  The  advantage  of  this  approach  is  the 
appearance  of  non-dimensional  parameters  such  as  Mach  number,  Reynolds  number, 
and  Prandtl  number,  which  can  then  be  varied  independently.  Denoting  dimensional 
variables  with  a  *,  the  nondimensionalization  appears  as  follows  (6:191)  (5:16) 
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Appendix  B.  Jacobian  Matrices 


The  Jacobian  matrices  employed  in  equations  (98-99)  are  presented  below. 
They  are  obtained  by  rewriting  the  vectors  E\,  Ei,  Vi,  and  Wi  appearing  in  equar 
tions  (51-56)  in  terms  of  the  conserved  variables  (pfJ,  pa/J,  pvfJ,  ptfJ )  and  dif¬ 
ferentiating  with  respect  to  U,  Uf,  or  as  indicated  in  equations  (98-99).  The 
resulting  matrices  are  then  (18:45) 
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where  <j>  =  1/2(7  —  1)(«3  +  v3)- 
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The  Jacobian  matrices  M  and  N  are  defined  by 
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The  coefficients  6,  and  d,  appearing  in  M  and  N  respectively  are  defined  in  equations 
(59)  through  (71). 
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Appendix  C.  Boundary  and  Initial  Conditions 


The  particular  solution  of  all  differential  equations  requires  the  specification  of 
boundary  and/or  initial  conditions.  In  this  case,  the  boundary  conditions  take  the 
form  of  specifying  flow  variables  on  the  airfoil  surface  and  along  the  grid  boundaries. 
Initial  conditions  are  specified  by  assigning  flow  variables  throughout  the  held  upon 
code  initiation. 

Refering  to  Figure  1,  freestream  values  in  non-dimensional  form  are  assigned 
throughout  the  flow  field  and  along  the  boundary  ABC. 

c  =  =  0 

Uoo  =  1  (119) 

v  ss  (Jo,  sin(a). 

The  initial  eddy  viscosity  is  generally  set  equal  to  a  constant  in  the  boundary  along 
the  airfoil.  Since  the  initial  conditions  are  a  function  of  Mach  number  and  angle 
of  attack,  it  is  required  that  restart  solutions  maintain  the  same  Mach  number  and 
angle  of  attack. 

At  the  airfoil  surface,  GFG,  an  adiabatic  condition  is  imposed.  The  velocities 
are  determined  there  by  no-slip  condition. 


P  =  Poo  =  1 
P=Poa= 

«  =s  Uoo  cos(a) 


«»  =  */*»  =  0. 


(120) 


The  bondary-layer  assumption  is  employed,  as  the  gradient  of  pressure  normal  to 
the  surface  is  assumed  zero.  The  adiabatic  assumption  implies  that  the  temperature 
gradient  normal  to  the  surface  is  zero.  The  surface  density  is  permitted  to  vary  as  a 
function  of  surface  pressure  and  temperature  via  the  perfect  gas  law. 
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PXM  —  PV(P,T). 


(121) 

(122) 


The  downstream  outflow  conditions  along  boundaries  AO  and  HC  are  com¬ 
puted  by  the  following  means 


P  =  Poo 


(123) 


±  ‘  =0. 
ae  « 


(124) 


In  general,  the  derivative  conditions  here  allow  for  large  scale  flow  fluctuations  due 
to  the  forces  applied  on  the  body. 

In  the  case  of  the  C-grid,  the  boundary  conditions  along  the  branch  cut  must  be 
specified  to  assure  continuity  is  satisfied.  Thus  along  cut  ED  and  EH,  the  following 
is  applied 


(125) 


102 


which  reduces  to  averaging  across  the  wake 
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(126) 


In  the  case  of  the  O-grid,  the  placement  of  the  wake  cut  is  arbitrary  and  in  this 
case  is  chosen  to  extend  from  the  leading  edge  to  the  forward  boundary.  Coincident 
grid  lines  are  superimposed  over  those  lines  bordering  the  O-grid  wake  cut  to  allow  for 
central  differencing  to  be  applied  across  the  wake  rather  than  the  somewhat  artificial 
boundary  condition  imposed  in  the  C-grid  wake.  This  represents  a  difference  in  the 
manner  in  which  the  boundary  conditions  are  applied  in  the  two  cases. 
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Appendix  D.  Turbulence  Model 


One  of  the  most  difficult  problems  in  modem  computational  fluid  dynamics  is 
the  correct  modeling  of  the  turbulence.  The  turbulence  model  employed  here  is  a 
modified  version  of  the  algebraic  eddy-viscosity  model  of  Baldwin  and  Lomu  (3). 
The  application  of  the  turbulence  model  is  applied  in  three  seperate  regions:  the 
boundary  layer  on  the  airfoil,  the  wake  proximal  to  the  airfoil  trailing  edge,  and  the 
far-field  wake  (see  Figure  59). 


The  boundary  layer  of  the  airfoil  is  divided  into  an  inner  and  outer  region.  In 
the  inner  region,  the  eddy  viscosity,  eit  is  given  by  the  Prandtl-Van  Driest  formula 
(18:12) 


with 


and 


e  =  p(ICD)7  |W| 

(127) 

"-"-[-i  ml 

(128) 

du  dv 

(129) 

W  dy  dx 

where  w  is  the  vorticity,  Y  represents  the  normal  distance  to  the  airfoil  surface, 
K  as  0.40  is  the  von  Karman  constant,  and  the  subscript  w  denotes  conditions 
located  at  the  surface  (wall). 

In  the  outer  region  of  the  boundary  layer,  the  eddy  viscosity  is  defined  by 
(18:12) 


e0 

Fk 


pKCcpYmn.Fm^F, 


(130) 

(131) 
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where  Fmax  =  max(Y|u;|D),  Ymas  is  the  value  of  Y  corresponding  to  Fm„  and  k  = 
0.0168,  =  1.6,  Ck  =0.3. 

The  turbulence  model  switches  from  the  inner  to  outer  formulation  at  the  first 
value  of  Y  for  which  c,  >  e0.  Transition  from  laminar  to  turbulent  flow  is  specified 
by  user  input  and  is  chosen  to  match  boundary-layer  trip  locations  when  comparing 
to  experimental  data. 

The  turbulence  model  in  the  far-wake  is  modeled  by  defining  the  turbulent 
eddy  viscosity  there  as  (18:13) 

Yma*AU2 

€vk  =  pCvk—jjr - Ft  (132) 

*  max 

where 

Cmk  ~  0.058 
Fmas  =  m«*(y|w|) 

A U  =  +  ^ 

In  the  wake,Ymal  is  measured  from  the  wake  center!  ie  as  determined  by  the 
location  of  minimum  velocity.  The  intermitancy  factor  Fk  is  obtained  from  equation 
(131).  The  constants  Cw±  and  C*  are  chosen  to  match  theoretical  values  for  an 
incompressible  turbulent  wake  (14). 

In  the  near-wake  the  turbulent  eddy  viscosity  is  determined  by  exponentially 
smoothing  the  eddy  viscosity  profile  at  the  trailing  edge  to  the  far-wake  profile.  The 
following  functional  form  is  employed  (18:14): 

enw  —  2  {[«(**«»  ^0  e»k(&o>  Y)) 

+A['{x0,Y)-'(su,Y))} 


(133) 


The  distance  x0  —  x1t  is  chosen  to  be  of  order  106,  where  6  is  the  average  boundary 
layer  thickness  at  the  trailing  edge.  It  should  be  noted  that  Visbal  considers  the 
near-wake  formulation  a  preliminary  approach  within  the  context  of  the  algebraic 
eddy  viscosity  model 


Figure  64.  Regions  of  Application  for  Turbulence  Model 


Appendix  E.  Computer  Resources  and  Codes 


The  Navier-Stokes  code  was  executed  on  two  computer  systems.  The  primary 
computer  utilized  was  the  Aeronautical  Systems  Division  Cray  XMP,  located  at 
Wright- Patterson  Air  Force  Base,  Ohio.  Selected  data  runs  were  accomplished  on  the 
Air  Force  Institute  of  Technology  Stellar  computer,  also  located  at  Wright-Patterson 
Air  Force  Base.  Table  3  compares  execution  times  of  the  two  computers  (CPU  time). 
The  times  quoted  are  for  500  iterations  using  a  119x40  C-Grid.  The  computations 
on  the  Stellar  were  performed  using  single-  precision  (32  bit)  accuracy,  while  the 
runs  performed  on  the  Cray  XMP  were  performed  using  double-precision  (64  bit) 
accuracy.  The  front-end  machine  used  to  interface  with  the  Cray  XMP  was  the 
Aeronautical  Systems  Division  Cray-Vax. 

Additional  computer  resources  employed  included  the  FDL  Iris  and  Prime 
computers.  The  Iris  workstation  was  utilized  for  grid  generation  and  data  reduction. 
Data  plotting  was  accomplished  using  the  Prime  computer,  in  conjunction  with  the 
Talaris  laser  printer. 

A  computer  listing  is  provided  for  the  metric  analysis  program  used  to  calculate 
the  Jacobian,  aspect  ratio,  orthogonality,  and  skewness  parameters  for  a  given  grid. 
The  transformation  metrics  are  extracted  directly  from  the  Navier-Stokes  code,  the 
above  parameters  are  then  calculated  from  the  metrics.  Listings  of  the  Navier-Stokes 


|  computer 

CPU  time  (sec) 

C-Grid 

Iterations 

Ci 

Ch 

|  Stellar 

60.1  -~1 

119x40 

500 

0.35398 

0.029575 

I  Cray  XMP 

64.4 

119x40 

500 

0.35416 

0.029480 

Table  3.  Comparison  of  CPU  times:  Cray  vs  Stellar 
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code  for  the  C-grid  and  the  associated  data  reduction  codes  may  be  found  by  refering 
to  Boyles  (5).  Both  the  C-  and  O-grid  Navier-Stokes  codes,  and  the  grid  generation 
codes,  are  archived  on  the  Stellar  computer. 
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